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£NJ ■ In this letter we present a high-precision evaluation of the expansions in e — (3 — d)/2 of (up to) four- 

CO ' loop scalar vacuum master integrals, using the method of difference equations developed by Laporta. 

We cover the complete set of fully massive master integrals. 
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^ I 1 Introduction 

Oh, Higher-order pcrturbative computations have become a necessity in many areas of theoretical physics, be 

it for high-precision tests of QED, QCD and the standard model, or for studying critical phenomena in 
condensed matter systems. 

Most recent investigations employ a highly automated approach, utilizing algorithms that can be im- 
plemented on computer algebra systems, in order to handle the growing numbers of diagrams as well as 
jJJ ■ integrals which occur at higher loop orders. 

Computations can be divided into four key steps. First, the complete set of diagrams including symmetry 
factors has to be generated. For a detailed description of an algorithm for this step for the case of vacuum 
topologies, see PQ. Second, after specifying the Feynman rules, the color- and Lorentz-algebra has to 
be worked out. Third, within dimensional rcgularization, massive use of the intcgration-by-parts (IBP) 
technique [2] to derive linear relations between different Feynman integrals in conjunction with an ordering 
prescription can be used to reduce the (typically large number of) integrals to a basis of (typically a few) 
master integrals [3]. Practical notes as well as a classification of vacuum master integrals is given in 
Fourth, the master integrals have to be solved, either fully analytically, or in an expansion around the 
space-time dimension d of interest. 

It is the fourth step that we wish to address here. While most work has been and is being devoted to 
d = 4, perturbative results in lower dimensions are needed for applications in condensed matter systems, 
as well as in the framework of dimensionally reduced effective field theories for thermal QCD, where recent 
efforts have made four- loop contributions an issue [5]- 

A very important subset of master integrals are fully massive vacuum (bubble) integrals, since they 
constitute a main building block in asymptotic expansions (see e.g. [3]). They are also useful for masslcss 
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theories, when a propagator mass is introduced as an intermediate infrared regulator \7\. 

The main purpose of this note is to numerically compute the complete set of fully massive vacuum master- 
integrals in terms of a high-precision e-expansion in d = 3 — 2e dimensions, in complete analogy with the 
four-dimensional work of S. Laporta [H]- 

The plan of the paper is as follows. In Section [3 we give a brief review of the method of difference 
equations applied to vacuum integrals. In Section^ we discuss the actual implementation of the algorithm. 
In Section 0] we display our numerical results for the truncated power series expansions in e of all fully 
massive master integrals, up to four-loop level, in d — 3 — 2e. 



2 The evaluation of master integrals through difference equations 

The method we have chosen to compute the coefficients of the truncated power series expansions of the mas- 
ter integrals is based on constructing difference equations for the integrals and then solving them numerically 
using factorial series. This approach was recently developed in Ref. [2], and below we briefly summarize its 
basic concepts following the notation of the original paper, which contains a much more detailed presenta- 
tion on the subject. While the method is completely general as it applies to arbitrary kinematics, masses 
and topologies 9 , our brief summary is somewhat adapted to the specific case of massive vacuum integrals. 
The main idea is to attach an arbitrary power x to one of the lines of a master integral U , 

where the Di = (pf + 1) denote inverse scalar propagators. In our case all of these share the same mass m, 
which we have therefore set to 1, noting that it can be restored in the end as a trivial dimensional prefactor 
of each integral. The original integral is then just U = U(l). Depending on the symmetry properties of 
the integral, there can be different choices for the 'special' line with the arbitrary power x, but in the limit 
x = 1 they all reduce to the original integral U. This degeneracy can (and will later) be used for non-trivial 
checks of the method. 

Employing IBP identities in a systematic way, it is possible to derive a linear difference equation obeyed 
by the generalized master integral U(x), 

R 

Y,Pi(*)U{x + 3) = F(x), (2) 
j=o 

where R is a finite positive integer and the coefficients pj are polynomials in x (and the space-time dimension 
d). The function F on the r.h.s. is a linear combination of functions analogous to U{x) but derived from 
simpler master integrals, i.e. integrals containing a smaller number of loops and/or propagators. 

The general solution of this kind of an equation is the sum of a special solution of the full equation, Uo(x), 
and all solutions of the homogeneous equation (F — 0), 

R 

U(x) = U (x) +J2u j (x), (3) 



where each (j — 0, R) 



°° V(r 4- 1 1 



is a factorial series 1 . Substituting into Eq. (J2J, one obtains the coefficients fi and K (the latter being a 
function of d), as well as recursion relations for the ^-independent coefficients a(s) (being functions of d as 

1 For a rigorous definition of the concept as well as a motivation for this kind of an ansatz, we refer the reader to Ref. 
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well) for each solution. For the homogeneous solutions, these recursion relations relate all coefficients to 
their value at s — 0, <ij(s) = Cj(s) a.j(0), where the Cj(s) are rational functions (of d as well). For the special 
solution, the ao(s) are completely fixed in terms of the inhomogeneous part F(x), consisting of 'simpler' 
integrals which are assumed to already be known in terms of their factorial series expansions. 

What remains to be done is to fix the x- and s-independent constants Oj(0), j ^ 0, in order to determine 
the weights of the different homogeneous solutions. To this end, it is most useful to study the behavior of 
U(x) at large x, where the first factor in 

U{x) = I RTTF 9{pi) (5) 

peaks strongly around p\ — 0. Hence, the large- a; behavior of the modified master integral is determined by 
the small- momentum expansion of the two-point function g(pi), which has one loop less than the original 
vacuum integral. In fact, for all cases we cover here, the first coefficient in the asymptotic expansion suffices. 
This is furthermore particularly simple, since it factorizes into a one-loop bubble carrying the large power 
x and a lower-loop vacuum bubble <?(0), which corresponds to U(x) with its 'special' line cut away, 

lim U (x) = 

x — >oo 

A comparison with the large-x behavior of Eqs. Q, (J3J, proportional to J^j ( J >j a j(fy x > can now be used 
to fix the Oj(0), of which maximally one will turn out to be non-zero for our set of integrals. 

Having the full solution at hand, we have in principle completed our entire task, as in the limit x — 1 
we recover from U (x) the value of the initial integral. Let us, however, add a couple of practical remarks 
here. What is still to be done is to perform the summation of the factorial series of Eq. which means 
truncating the infinite sum at some s max . Studying the convergence behavior of these sums, one notices 
that even in the cases where they do converge down to i ~ 1, their convergence properties usually strongly 
decline with decreasing x. This means that in practical computations, where one aims at obtaining a 
maximal number of correct digits for U(l) with as little CPU time as possible, the optimal strategy is to 
evaluate the integral U (x) with the factorial series approach at some :r max 1 and then use the recurrence 
relation of Eq. J5J to obtain the desired result at x — 1. The price to pay is, however, a loss of numerical 
accuracy at each 'pushdown' (x — » x — I ) step due to possible cancellations, which makes the use of a very 
high x max impossible. In practice the strategy is to determine an optimal value for the ratio s max /x max . 
To give an example, for the four-loop integrals of Section^ we have found that s max /x max ~ 50 is a good 
value, while we used a range of s max ~ 1350 . . . 2000. 

3 Implementation of the algorithm 

As is apparent from the preceding section, there are three main steps involved in obtaining the desired 
numerical coefficients in the e-expansion of each master integral: deriving the difference equations obeyed 
by each integral, solving them in terms of factorial series, and finally performing the e-expansion and 
numerically evaluating the sum of Eq. (0J (truncated at s max ) to the precision needed. We will briefly 
address each of them in the following. 

For the first step, we slightly generalized the IBP algorithm we had used for reducing generic 4-loop bubble 
integrals to master integrals, which follows the setup given in and whose implementation in FORM ^D] 
is documented in 0|. The main difference is an enlarged representation for the integrals, keeping track of 
the line which carries the extra powers x, as well as the fact that there are now two independent variables 
(d, x), requiring factorization (and inversion) of bivariate polynomials, as opposed to univariate polynomials 
in the original version. 

Second, staying within FORM for convenience, we implemented routines that straightforwardly solve 
the difference equations in terms of factorial series, along the lines of |3J- This is done starting with the 
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simplest one-loop master integral, and working the way up to the most complicated (most lines) four- 
loop integral, ensuring that at each step, the 'simpler' terms constituting the inhomogeneous parts of the 
difference equation are already known. The output are then plain ascii files specifying each solution in the 
form of Eq. (J3J as well as containing recursion relations for the coefficients a(s). Note that these first two 
steps are performed exactly, in d dimensions. 

Third, once the recursion relations for the coefficients a(s) were known, we used a Mathematica program 
to obtain their numerical values at each s to a predefined precision, and to perform the summation of the 
factorial series. While this procedure is in principle very straightforward, there are some twists that we 
employed to help reduce the running times significantly, most of which are probably quite specific to our 
use of Mathematica. To avoid a rapid loss of significant digits in solving the recursion steps that relate each 
a(s) to a(0), especially those for the homogeneous coefficients, we first solved the relations analytically and 
only in the end substituted the numerical value (actually the truncated e-expansion) of the first non-zero 
coefficient. In fact, we found Mathematica to operate quite efficiently with operations like multiplication of 
two truncated power series, so that we relied heavily on it. Furthermore, since — not surprisingly — the 
most time-consuming part in the summation of the series turned out to be the e-expansion of T-functions, 
we achieved a notable speed-up by substituting the T-functions with large arguments by suitable products 
of linear factors times T- functions of smaller arguments. Finally, a vital step in avoiding an excessive loss in 
the depth of the e-expansions when going from one integral to the next, was to apply the 'Chop' command 
to remove from the results and coefficients excess unphysical poles, whose coefficients were of the order of, 
say, 10 -50 or less. 

4 Numerical results 

Below we list the Laurent expansions in e = (3 — d)/2 of the 1+1+3+13 fully massive vacuum master 
integrals up to four loops. We use an intuitive graphical notation, in which each line represents a massive 
scalar propagator, while dot on a line means it carries an extra power. The integral measure we have chosen 
here is 



1 /-d 3 - 2e 



?■ (7) 



r(3/2 + e) J it 3 / 2 -" ■ 

In each case 2 we provide the first 8 e-orders keeping the accuracy at 50 significant digits for the 1-, 2-, and 
3-loop master integrals and at 22-25 for the 4-loop ones. To obtain more e-orders and significant digits is 
merely a matter of additional CPU time. 



o 



e 



4.0000000000000000000000000000000000000000000000000 

16.000000000000000000000000000000000000000000000000 e 2 
64.000000000000000000000000000000000000000000000000 e 4 

256.00000000000000000000000000000000000000000000000 e 6 + O (e 8 ) (8) 

4.0000000000000000000000000000000000000000000000000 e" 1 
14.487441729730630111648209847429586185151846775400 
41.495035953369978394225958244504121655360756728405 e 
107.49752321579967383991953818365893067117808339742 e 2 



_> 



With the exception of the last two integrals, for which we were at this time able to produce only the first 6 and 5 e-orders, 



respectively. 
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263.49878761720606330238135348797499506915058750280 e 3 
623.49940078392000186832902635721463645559035022216 e 4 
1439.4997026869879573968449524699557874962297882621 e 5 

3263.4998520860644726225542919399943895943491031166 e 6 + O (e 7 ) (9) 

64.00000000000000000000000000000000000000000000000 e" 1 
49.44567822334599921081 142309329320142732803439623 
1981.207736229513534030093683214422278348416661525 e 
235.7077170926718752095474374908098006136204356228 e 2 
63521.71508871044639640714223746384514019533126715 e 3 
33675. 1 1 1 1 1780076696716334804652776927940758434016 e 4 
2213147.071275511251113640247844877948334091419700 e 5 

1414250.728717593474053272387541196652013773984236 e 6 + O (e 7 ) (10) 

32.859770043923503738827172731532536947448547448996 
365.41238154175547388711920818936800707879030719734 e 
2803.7940402523167047150293858439985472095966118207 e 2 
18727.187392108144301607279844058527418378836943988 e 3 
117794.35873133306139734878960626307962150043480498 e 4 
721386.63300305569920915438185951112611780543107044 e 5 
4366100.1639736899128559563097848872427318803864139 e 6 
26291285.708454833832306242766439811661977583440814 e 7 + O (e 8 ) (11) 

1.391204885296021941812048136925327740910466706390 
4.898152455251800666032641168608190942446944333758 e 
12.98842503803858164353982398007130232261458098462 e 2 
30.39637625288207454078370310227949470365033235457 e 3 
66.67957617359017942652215661267829752624475575093 e 4 
140.9974945708845413812214824315460314748605690042 e 5 
291.7287632268179138442199742398614147733926624689 e 6 

595.7006275449402266695675282375932229509102799733 e 7 + O (e 8 ) (12) 

720.0000000000000000000000 c" 1 - 52.13034199729620858728708 
33748.69042965137616701638 e + 10819.60558535024688749473 e 2 
1311729.690542895866693548 e 3 + 615270.7589383441011319577 e 4 

48899219.67276170476701364 e 5 + 24885879.11003549349511900 e 6 + O (e 7 ) (13) 

32.00000000000000000000000 c" 1 + 21.28521367989184834349148 
945.4764617862257950102533 e - 500.9879407913869195081538 e 2 
29027.99548541518650323471 e 3 - 34796.65982174097113175672 e 4 
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- 993306.5068744076465770453 e 5 - 1406349.173668893367086333 e e + O (e 7 ) (14) 

+ 8.000000000000000000000000 <T 2 - 25.94976691892252044659284 e" 1 

- 152.5193565764658289654545+ 2653.873458838396323815566 e 

- 23471.05910309626447406639 e 2 + 169839.2007120049515774452 e 3 

- 1124117.877397355450165203 e 4 + 7116455.837989754857686241 e 5 + O (e 6 ) (15) 

+ 78.95683520871486895067593 <T X - 1062.608419332108844057560 
+ 9340.076804859596283223881 e - 68699.47293187699594375521 e 2 
+ 462145.6926820632806821051 e 3 - 2963063.672524354359852913 e 4 

+ 18494675.22629230338091457 e 5 - 113673206.9834859509114931 e e + O (e 7 ) (16) 

+ 33.05150971425671642138224 - 358.4595946559340238066389 e 
+ 2451.469078369636793421997 e 2 - 13564.14170819716549262162 e 3 
+ 66602.55178881628657891800 e 4 - 303915.1384697444382333780 e 5 

+ 1323370.670112542076081095 e 6 - 5589978.086026239748023404 e 7 + O (e 8 ) (17) 

+ 27.57584879577521927818358- 291.4075344540614879796315 e 

+ 1956.162997112043390446958 e 2 - 10678.5639091187201818981 e 3 
+ 51925.3888799007705970928 e 4 - 235296.36309585614167636 e 5 

+ 1019555.9650538012793966 e 6 - 4292011.3101269758990557 e 7 + O (e 8 ) (18) 

+ 19.84953756526739935782082 - 200.9768306606422068619864 e 

+ 1308.883448000100198800887 e 2 - 6990.22562100063537185149 e 3 
+ 33456.8326902483214417013 e 4 - 149903.697032731221510018 e 5 

+ 644404.61801211590204150 e 6 - 2697912.0878890801856234 e 7 + O (e 8 ) (19) 

+ 3.141336279450209755917806 - 19.78740273338730374386071 e 

+ 83.81604328128850410126511 e 2 - 295.3496021971085625102731 e 3 
+ 934.2247995435558122394582 e 4 - 2751.31852347627462886909 e 5 

+ 7700.18972963585089750348 e 6 - 20740.9769474365145116212 e 7 + O (e 8 ) (20) 

+ 2.012584635078182771827701 - 10.76814227797251921324485 e 

+ 39.40636857271936487899035 e 2 - 121.0015646826735646109733 e 3 
+ 335.6942965583773421544251 e 4 - 872.009773755552224781319 e 5 

+ 2163.88707221986880315576 e 6 - 5193.51249188593850483093 e 7 + O (e 8 ) (21) 

+ 1.27227054184989419939788 - 5.67991293994853579036683 e 

+ 17.6797238948173732343788 e 2 - 46.5721846649543261864019 e 3 
+ 111.658522176214385363568 e 4 - 252.46396390100217743236 e 5 
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+ 549.30166596161426941705 e 6 - 1164.5120588971521623546 e 7 + O (e 8 ) (22) 

= +0.297790726683752651865168- 0.709896385699143430126726 e 

+ 1.40535549472683132370135 e 2 - 2.45721908509256673440117 e 3 
+ 4.00998036005764459707090 e 4 - 6.2518071963546459390185 e 5 

+ 9.4402506572040685160665 e 6 - 13.924465979877416801887 e 7 + O (e 8 ) (23) 

= + 0.233923932580303206470057 - 0.48523164074102176840584 e 
+ 0.88555744401503729577888 e 2 - 1.438019871368410241810 e 3 
+ 2.198725350440790755608 e 4 - 3.231974794381719679729 e 5 + O (e 6 ) (24) 

= + 0.195906401341238799905792 - 0.37006152907989745845214 e 
+ 0.65228273818146302130509 e 2 - 1.029288152514143871118 e 3 
+ 1.542484509438506710808 e 4 + O (e 5 ) (25) 

We have performed various checks in order to test the correctness of our recursion relations as well as to 
verify the number of exact digits contained in our results Eqs. JSJ-lj^SJ. The first task we have completed 
by exploiting the fact that the recursion relations are not specific to d = 3 — 2e, but can easily be applied to 
any dimension, such as d = 4 — 2e. We have successfully verified the results of Ref. 5] to somewhat lower 
accuracy and depth in e. Note that our choice of a basis for 4-loop master integrals differs slightly from the 
one made in [Sj. The relations needed for a basis transformation are listed in |3J. An immediate advantage 
in the light of difference equations is that with our choice, the above results Eqs. 114ll . (|20|l and (|23|l follow 
'for free' from their counterparts without dots. 

The accuracy of our three-dimensional results we have on the other hand examined in three independent 
ways: 

• by comparing the numerical results to existing analytic calculations; they can be found in |11| (diver- 
gent and constant parts of Eqs. ©-(EJ), H2 (leading term of Eq. [EH (divergence of Eq. lfTB)ll 
and |14I15| (all divergences and some constant parts of 4-loop integrals, as well as some 0(e) terms of 
lower-loop cases). 

• by comparing the results obtained by raising topologically incquivalcnt lines to the power x, 

• by analyzing the convergence properties of the factorial series, i.e. by checking the stability of our 
results with respect to varying s max . 

The first method is of course exact, but is only available for a few low (in e) orders for approximately 
half of the integrals considered. The second one, on the other hand, has the advantage of covering all the 
different powers of e, but is inapplicable for those integrals, in which all propagators are equivalent (e.g. the 
basketball-topology). The third method is then the most widely applicable one, but has the downside of 
providing no evidence for the correctness of our results, rather giving only the number of digits stable in the 
variation of the cut-off of the factorial series. For the integral of Eq. (|25|l only the last method is available, 
but in addition we have verified the leading term in the result to 3 digits using a Monte Carlo integration 
of an 8-dimcnsional integral representation derived for this integral in Ref. |14| . 

One might be concerned about the rapid growth with increasing e-orders of most of the coefficients. This 
is, as was pointed out in [5], caused by poles that the integrals (seen as functions of d) develop near d = 3, 
e.g. at d = 7/2, 4, etc. It is to be expected that factoring out the first few of these nearby poles in each case 
will improve the apparent convergence in e considerably. 
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In principle, having a method at hand that is capable of generating coefficients to very high accuracy, 
even to a couple of hundred digits, one could now use the algorithm PSLQ 16 combined with an educated 
guess of the number content of some of the yet-unknown constant terms, in order to search for analytic 
representations of the numerical results. These could then in turn be used as an inspiration to find useful 
transformations of the integral representation of the original integral, which might allow for a fully analytic 
solution in those cases where it could not yet be achieved. We have not made any attempts in that direction, 
since the numerical accuracy of the results Eqs. ^-1)25(1 should be sufficient for all practical purposes. 
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